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Abstract 

In this work, we study the critical behavior of an epidemic propagation model that considers 
individuals that can develop drug resistance. In our lattice model, each site can be found in one 
of four states: empty, healthy, normally infected (not drug resistant) and strain infected (drug 
resistant) states. The most relevant parameters in our model are related to the mortality, cure 
and mutation rates. This model presents two distinct stationary active phases: a phase with co- 
existing normal and drug resistant infected individuals and an intermediate active phase with only 
drug resistant individuals. We employ a finite-size scaling analysis to compute the critical points 
the critical exponents ratio (3/u governing the phase-transitions between these active states and 
the absorbing inactive state. Our results are consistent with the hypothesis that these transitions 
belong to the directed percolation universality class. 

PACS numbers: 64.60.Ht, 05.70.Ln, 87.23. Cc, 87.19.Xx 
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I. INTRODUCTION 



The spread of epidemic processes in a population has been profusely studied by biolog- 
ically oriented mathematicians mainly by the use of a set of partial differential equations 
governing the time evolution of global variables such as the density of infected, refractory 
and healthy individuals At high dimensions where local density fluctuations are neg- 
ligible, this mean-field like approach gives the correct overall behavior of these systems, 
according to renormalization group arguments. At low dimensions, density fluctuations be- 
come relevant and microscopic approaches have to be considered. In general, most models of 
epidemic propagation processes present a second order phase transition between a vacuum 
and a stationary active state 0,0,0,0,0]. The vacuum corresponds to an absorbing state 
towards which the system can evolve but from which it cannot leave. Systems with such 
absorbing states cannot satisfy detailed balance due to the irreversible character of their 
dynamics. 

Nonequilibrium phase transitions from active fluctuating into inactive absorbing states 



have been studied extensively for the last twenty years |3l|. These phase 

transitions into absorbing states are characterized by a nonequilibrium critical behavior 
similar to that of an equilibrium phase transition in many respects. One usually uses the 
concept of scale invariance to understand nonequilibrium phase transitions with critical 
exponents characterizing its universal behavior. These exponents allow us to categorize 
nonequilibrium phase transitions into different universality classes. In reaction models for 
the spread of epidemics, the density of infected individuals pi acts as the order parameter, 
vanishing in the absorbing state and having a finite average value in the active state. The 
transition between the active and stationary state is achieved by tuning a typical reaction 
rate A that acts as the relevant control parameter. In the vicinity of the critical point, the 
density of infected individuals vanishes as p oc (A — \c)^ . Nonequilibrium phase transitions 
are often characterized by spatial and temporal correlation lengths. Close to the transition 
point, the spatial and temporal correlation lengths diverge as oc (A — \c)~'^ and oc (A — 
Ac)^*^'. These two correlation length scales are related by = where z = u/u' is called the 
dynamic exponent. These three critical exponents ( f3, u, z) form the basic set of exponents 
characterizing the universality class of a given reaction-diffusion process. It is generally 
believed that nonequilibrium phase transitions into absorbing states can be categorized into 
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a finite number of universality classes. Usually, systems presenting a continuous transition 
to a single absorbing state belong to the directed percolation universality class ^, 11, 1^ . 
However, some exceptions to this general rule have been receritly reported to occur basically 



in systems with a infinite number of absorbing states 



diffusion 
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or with particle 



20l|21|]. 



In recent years, there has been a growing interest in studying diseases caused by organisms 
that can mutate into drug resistant forms. One of the major reasons for the emergence of 
such resistant organisms is the inappropriate use of the drugs used to combat them. Such is 
the case, for example, when antibiotics are incorrectly administered to treat the tuberculosis 



bacillus 
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or the excessive use of anti-malaria drugs in Amazonia, Brazil 
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28l |. Although these drugs are the only known effective treatments for such diseases. 



their abuse can very quickly lead to the emergence of drug resistant mutant strains. For 
example, while the exact mechanisms are not fully understood, strains have evolved that 
are resistant to every drug used in Amazonia during the 1950's and 60's. To retard the 
emergence of such resistant strains, measures are now being taken in the medical profession to 
administer drugs only when needed, and giving preference to drug "cocktails" with multiple 



components 
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28|. Recent 



mutations has been introduced 



y, an epidemic spreading model with immunization and 
29|. The model also exhibits a phase transition belonging to 
the directed percolation universality class and shows explicitly that the protection gained 
by immunization is substantially decreased by the occurrence of mutations. 

In this work, we will study a 1-dimensional stochastic epidemic propagation model consid- 
ering the potential emergence of a population of drug resistant organisms. We will be mainly 
interested in determining the stable phases and the scaling exponents governing the possible 
phase transitions. The model presents two coupled order parameters, namely, the density 
of drug sensitive and drug insensitive individuals. The coupling of an order parameter to 
a fluctuating density has been shown to modify the universalit v class of the absorbing sate 
phase transition in some models of population dvnamicsjsfll. |^. By employing extensive 
numerical simulations in finite populations and finite size scaling analysis, we will show that 
an active phase emerges in which there is no drug sensitive individuals but a finite density 
of drug resistant ones. The relative strength of the parameters controlling the processes 
of healing, death and development of drug resistance, actually determines the stationary 
population state. Further, we will show that all identified transitions in this model seem to 
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belong to the directed percolation universality class. 



II. BASIC MODEL FOR THE EMERGENCE OF DRUG RESISTANCE 

Here, we will follow closely the ideas of [s^. Individuals distributed at random on a 
given environment can be in three possible states, namely healthy, infected by normal strain 
and infected by resistant strain. Birth of individuals will be introduced by considering that 
a vacant site becomes occupied by a healthy individual at a rate j3i. Infection occurs by 
contact at a rate /SiUi where Ui is the number of infected neighbors and the index i stands 
for normal (z = 2) and drug resistant {i = 3) states. Individuals infected with the normal 
strain can heal at a rate pc- In what follows, we will consider that the healing rate is due 
to drug administration, i.e., the spontaneous healing rate is negligible. Infected individuals 
can also develop drug resistance at a rate pr- This process can be originated by both natural 
mutations and bad drug administration. Individuals that develop drug resistance will not 
be allowed to heal and die at a mortality rate pm- The average surviving-time of drug 
resistant individuals is much smaller than the average life-time of individuals infected with 
the normal strain, which will be taken as infinite in our simulations. The above dynamical 
rules describing the processes of birth, death, healing, development of drug resistance and 
infection by contact with normal or resistant strains are summarized in table 1. Notice that 
the state with only healthy individuals is an absorbing one, since it is no longer possible for 
anybody to become infected. Such is the case, for example, with small pox, which has been 
eradicated. 

A mean field approximation can be described using the following rate equations for the 
average densities: 

dpi 
dt 
dp2 



= PlPO - P2P1P2 - PsPlPS + Pcp2 



P2P1P2 - Prp2 - Pcp2 
^ = /^sPlPS + Prp2 - Pmp3 

at 

1 = 



where a homogeneous, i.e., spatially uniform density of all populations was assumed. It is 
straightforward to show that there are three possible stationary states: (i) the absorbing 



state with only healthy individuals; (ii) an active state where healthy individuals coexist 
with individuals infected by normal and resistant organisms; and (iii) a second active state 
where healthy and infected individuals coexist but all infected individuals have the resis- 
tant organism. An important point to note is that normal organisms always coexists with 
resistant organisms, i.e., a stationary active state without a drug resistant population is not 
possible. 

The basic contact process (CP), introduced by Harris js^ as a "toy model" for sprea ding 
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of disease, corresponds to the limiting case where only states 1 and 2 are allowed 
0]. The infection reaction rate is An where n is the number of infected neighbors, and the 
cure rate is 1, i.e. all sick individuals will recover in the next time step. Below the critical 
threshold A < Ac, the absorbing state is reached, but for A > Ac a stationary active state 
is possible. The model can also be applied to study the spread of a disease in the absence 
of mutations. In this case, we have only states 0, 1 and 2. This model is known in the 
literature as the forest fire model The above two limits have been extensively 

investigated in the literature and accurate estimates of their critical exponents have shown 
that they belong to the directed percolation universality class. Here, we are going to perform 
extensive numerical simulations to obtain the full phase-diagram of the model allowing all 
possible states. 



III. NUMERICAL SIMULATIONS: PHASE DIAGRAMS AND CRITICAL EX- 
PONENTS 

We employed the above described model in finite one-dimensional chains with periodic 
boundary conditions. Sites are considered to be vacant or occupied by an individual that 
interacts with its nearest neighbors. The initial state corresponds to a random distribution 
of healthy {rji) and infected (172) individuals. For a given set of model parameters, we let 
the system evolve towards a statistically stationary state. In what follows we are going to 
restrict our study to the particular case of /5i = = /^a = 1, letting as free parameters the 
rates of mutation (pr), cure (pc) and mortality (pm)- In our simulations, the statistically 
stationary state was considered as reached after = 10^ lattice sweeps. Configurational 
averages were performed by considering the subsequent 2 x 10^ configurations. Finite size 
simulations of dynamical systems exhibiting an absorbing state are usually compromised 
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by the fact that this is the only stationary state. To avoid this feature, we reactivate 
the dynamics whenever the absorbing state is reached by introducing a small density of 
the critical order parameter. The dynamics follows the rules of contamination, healing, 
mutation, death and birth processes summarized in table 1. In addition to these rules, the 
contamination of a healthy individual that is surrounded by both kinds of infected sites is 
chosen to be preferentially updated to state 7^3 (resistant strain). 

In figure 1, we show the equilibrium density of normally infected individuals (^2) and 
drug resistant ones (ps) as obtained from simulations on chains with L — 800 sites for three 
characteristic values of the mutation rate pr- We observe that the absorbing state is the 
equilibrium state for large values of cure and mortality rates. Two distinct active stationary 
regimes are also clearly identified. The first one corresponds to a state with zero density 
of normally infected but with a finite density of drug resistant ones. This phase is stable 
at low death rates and its emergence is independent of the mutation rate 7^ 0. The 
large surviving time of drug resistant individuals in this regime favors the contamination 
process. The second active state has co-existing densities of normal and drug resistant 
infected individuals. It emerges in the regime of low cure rates but disappears for large 
mutation rates. A direct transition between these active states takes place at low death and 
cure rates. 

In order to closely investigate the transition between the active state with co-existing 
densities and the absorbing state, we computed the order parameter densities for pc — 
(no heahng), which is the most favorable case for co-existence. In figure 2, we show the 
relevant densities as function of pr and pm- Notice that there is a critical mutation rate 
above which the phase with co-existing densities disappears. The maximum value for the 
critical mutation rate is attained for = 1. To have an accurate estimate of the critical 
parameters, we employed a finite size scaling analysis of our numerical data. Figure 3 shows 
the density p2 as a function of the mutation rate for the special case of Pc = and Pm — i 
and several chain sizes. The data suggest a continuous transition near pr — 0.175. A 
numerical renormalization group analysis of the above data can provide a better estimate 
of the critical value. This procedure explores the finite size scaling relation that the order 
parameter density is expected to obey at the vicinity of the critical point, namely 

P2{L,Pr)=L-^/'^f[{Pr-p;)L'/''] . (1) 
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Therefore, one can construct a family of auxiliary functions 

According to the finite size scaling hypothesis, these auxiliary functions shall all intercept at 
the critical point and gip*) — —(i/^- In figure 4, we report the above scaling analysis to the 
data shown in figure 3. All curves obtained from distinct pairs of lattice sizes intercept nearly 
at the same point. Without considering any correction to scaling, we obtain an accurate 
estimate of p* = 0.176(1) and the critical exponents ratio 13/v = 0.25(1). The above 
finite-size scaling analysis was employed to obtain the full phase diagram in the absence of 
treatment (pc = 0) as shown in figure 5. The absorbing healthy state (H) is always the 
equilibrium one at large mutation and mortality rates. The phase diagram exhibits three 
distinct phase transitions. For mortality rates above p*^ — 0.224(1), the active state is 
the one with co-existing densities of normal (N) and resistant (R) individuals which is the 
stationary state at small Pr- For mortality rates below p^ the stationary state is always 
active, as the long surviving-time of infected individuals makes the contamination process 
the predominant one. However, there is a transition between the active co-existing state 
at small Pr to the active state with a finite density of drug resistant individuals and a 
zero density of normal individuals at large Pr- Notice that at intermediate mutation rates Pr 
the phase-diagram presents re-entrant transitions with the sequence of active-inactive-active 
equilibrium states appearing with increasing mortality rates Pm- In all cases, the measured 
values oi P/u were consistent with these transitions belonging to the directed percolation 
universality class. In figure 6 we illustrate the main phase-diagram for a fixed mutation rate. 
In this case, the absorbing state is the equilibrium state at large cure and mortality rates. At 
low mortahty rates [below the critical value = 0.224(1)], there is always a finite density 
of drug resistant individuals. With increasing cure rate Pc a transition between the active 
co-existing state and the absorbing state can takes place. The line delimiting this transition 
saturates at a maximum critical cure rate p*(pr = 0.05) = 0.219(1) which decreases with 
increasing mutation rates. Further, this co-existing phase is only stable for mortality rates 
larger than p5^(pr = 0.05) = 0.168(1) and moves to larger values as Pr increases. At small 
mutation rates the above two lines intercept (as shown in the figure) and a direct transition 
between the two possible active states can take place. On the other hand, at the limiting 
value of p* — 0.176(1) the NR to H transition fine disappears and the co-existing phase is 



no longer stable. Re-entrant sequence of active-inactive-active equilibrium states can also 
be achieved at intermediate cure rates. 

IV. SUMMARY AND CONCLUSIONS 

We have studied a stochastic epidemic model in a 1-dimensional lattice which considers 
the possibility of individuals to develop drug resistance. Individuals with the normal disease 
may heal due to drug administration but those with the resistant strain do not heal and 
die at a finite mortality rate. We numerically simulated this model of population dynamics 
in finite chains with periodic boundary conditions and investigated the properties of the 
stationary states. The model presents one absorbing inactive state with only healthy indi- 
viduals. There is also a main active stationary state with fluctuating densities of normal 
and drug resistant infected individuals. Further, an intermediate active stationary state 
with no normally infected individuals and a fluctuating density of drug resistant ones is also 
realizable. This state can be the stationary one only when the mortality rate is below a 
critical threshold — 0.224, since under this condition the long surviving-time of drug 
resistant individuals potentialize this epidemic process. For mortality rates above the crit- 
ical threshold a transition between the main active to the absorbing state takes place with 
increasing cure rates. We reported typical phase- diagrams showing the stationary states as 
a function of the main transition rates. The phase-diagram exhibits re-entrant regions and 
a direct transition between the main and intermediate active states at low mutation rates. 
We further employed a finite-size scaling analysis to compute the exponent ratio P/u which 
is consistent with these transitions belonging to the directed percolation universality class. 
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TABLE I: Summary of the relevant processes, their effects on the population and transition rates 
for the present epidemic model with drug resistance. 
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FIG. 1: The stationary densities of normally infected (^2) and drug resistant (ps) individuals as a 
function of cure rate (pc) and mortality rate (pm) for typical mutation rates. The three possible 
stationary states are clearly identified, namely, the absorbing state {p2 = Ps = 0), the main active 
state (finite p2 and ps) and the intermediate active state (p2 = and finite ps). The intermediate 
active state is practically independent of the cure and mutation rates. The main active state is 
only stable for small mutation rates and disappears for Pr > Pr = 0.176. 
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Pc = 0.0 




FIG. 2: The stationary densities of normally infected {p2) and drug resistant (ps) individuals as 
a function of mutation and mortality rates for the special case of Pc = (no healing). Here the 
transition from the main active to the absorbing state is better illustrated. In this phase a small 
density of drug resistant individuals co-exists with a large density of normally infected ones. 
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FIG. 3: The stationary density of normally infected individuals (^2) versus the mutation rate Pr 
for the limiting case of Pc = and Pm = 1. Results from simulations on chains with different sizes 
are shown. The continuous phase transition from the active to the absorbing state is smoothed by 
finite size effects. 
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FIG. 4: Numerical renormalization group of the same data from figure 3 in the vicinity of the 
transition point. Renormahzation using several pairs of chain sizes are shown. All auxiliary 
functions intercept at a common point which determines the critical value p* = 0.176(1) and the 
critical exponents ratio (3/iy = 0.25(1). This exponent is consistent with this transition belonging 
to the directed percolation universality class. 
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FIG. 5: Phase-diagram for the special case of Pc = 0. The absorbing state (H) is only stable for 
mortality rates above = 0.224(1). In this region, a transition from the main active (NR) to the 
absorbing state (H) occurs for increasing mutation rates. The absorbing state is reached trough the 
death process. Below pj^ there is a transition between the main active (NR) to the intermediate 
active state (R). 



16 
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FIG. 6: Phase-diagram for a typical mutation rate Pr = 0.05. For low mortality rates the stationary 
state is again always active and mostly with drug resistant individuals. The main active state may 
be stable at small cure rates. Above the critical mortality rate the transition from the active (NR) 
to the absorbing (H) state takes place with increasing cure rates. In this region the absorbing state 
is reached through both death and cure process. 
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